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Abstract 

Symplectic  Runge-Kutta  schemes  for  the  integration  of  general  Hamiltonian  systems  are  implicit.  In  practice  one  has  to 
solve  the  implicit  algebraic  equations  using  some  iterative  approximation  method,  in  which  case  the  resulting  integration 
scheme  is  no  longer  symplectic.  In  this  paper  we  first  analyze  the  preservation  of  the  symplectic  structure  under  two  popular 
approximation  schemes,  hxed-point  iteration  and  Newton’s  method,  respectively.  Error  bounds  for  the  symplectic  structure  are 
established  when  N  fixed-point  iterations  or  N  iterations  of  Newton’s  method  are  used.  The  implications  of  these  results  for 
the  implementation  of  symplectic  methods  are  discussed  and  then  explored  through  extensive  numerical  examples.  Numerical 
comparisons  with  non-symplectic  Runge-Kutta  methods  and  pseudo-symplectic  methods  are  also  presented. 


AMS:  65L05;  65L06;  65P10 

Key  words:  Geometric  integrators;  Hamiltonian  structure;  Symplectic  Runge-Kutta  methods;  Pseudo-symplecticness; 
Fixed-point  iteration;  Newton’s  method;  Convergence 


1  Introduction. 

Geometric  integration  methods  —  numerical  methods  that  preserve  geometric  properties  of  the  flow  of  a  differential 
equation  —  outperform  off-the-shelf  schemes  (e.g.,  fourth  order  explicit  Runge-Kutta  method)  in  predicting  the  long¬ 
term  qualitative  behaviors  of  the  original  system  (Hairer  et  ah,  2002).  For  systems  evolving  on  differentiable  manifolds 
(including  the  important  setting  of  Lie  groups),  geometric  integrators  that  preserve  the  manifolds  are  currently  a 
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subject  of  great  interest  to  theorists  and  practitioners.  See  for  instance  Budd  and  Iserles  (1999).  Applications  of  such 
techniques  are  of  interest  in  a  variety  of  physical  settings.  See  for  instance  Krishnaprasad  and  Tan  (2001)  for  results 
related  to  the  integration  of  Landau-Lifshitz-Gilbert  equation  of  micromagnetics. 

An  important  class  of  geometric  integrators  are  symplectic  integration  methods  for  Hamiltonian  systems.  See  Sanz- 
Serna  and  Calvo  (1994);  Marsden  and  West  (2001)  and  references  therein.  When  the  Hamiltonian  has  a  separable 
structure,  i.e.,  H{p,q)  =  T{p)  +  V{q),  explicit  Runge-Kutta  type  algorithms  exist  which  preserve  the  symplectic 
structure  (Forest  and  Ruth,  1990;  Yoshida,  1990;  Candy  and  Rozmus,  1991;  McLachlan  and  Atela,  1992).  However, 
for  general  Hamiltonian  systems,  the  symplectic  Runge-Kutta  schemes  are  implicit  (Sanz-Serna,  1988).  In  practice 
one  has  to  solve  the  implicit  algebraic  equations  for  the  intermediate  stage  values  using  some  iterative  approximation 
method  such  as  fixed-point  iteration  or  Newton’s  method. 

In  general,  with  an  approximation  based  on  a  finite  number  of  iterations,  the  resulting  integration  scheme  is  no  longer 
symplectic.  Error  analysis  on  the  structural  conservation,  like  the  analysis  on  the  numerical  accuracy,  provides  insight 
into  a  numerical  method  and  helps  in  making  judicious  choices  of  integration  schemes.  An  example  of  this  is  Austin 
et  al.  (1993),  where  the  error  estimate  for  the  Lie-Poisson  structure  was  given  for  integration  of  Lie-Poisson  systems 
using  the  mid-point  rule.  The  first  objective  of  this  paper  is  to  investigate  the  loss  of  symplectic  structure  due 
to  the  approximation  in  solving  the  implicit  algebraic  equations.  The  fixed-point  iteration-based  approximation 
and  Newton’s  method-based  approximation  are  analyzed,  respectively.  For  either  method,  an  error  bound  on  the 
symplecticness  of  the  numerical  flow  is  established  when  N  iterations  are  adopted  for  any  N  >  1.  It  turns  out  that, 
under  suitable  conditions,  the  convergence  rate  of  the  symplectic  structure  is  closely  related  (but  not  equal)  to  the 
rate  of  convergence  to  the  true  solution  of  the  implicit  equations.  Hence  the  methods  become  almost  symplectic  as 
N  gets  large. 

The  implications  of  the  error  bounds  for  implementing  symplectic  Runge-Kutta  schemes  are  then  studied  in  combina¬ 
tion  with  a  series  of  numerical  examples.  The  question  is  how  to  strike  the  right  balance  between  the  computational 
cost  and  the  structural  preservation.  Choice  of  the  step  size,  the  initial  iteration  value,  and  fixed  point  iteration 
versus  Newton’s  method  are  discussed.  Numerical  comparisons  are  also  conducted  with  non-symplectic  explicit 
Runge-Kutta  methods  and  with  pseudo-symplectic  methods  proposed  in  Aubry  and  Chartier  (1998).  Note  that 
pseudo-symplectic  integrators  are  explicit  and  designed  to  conserve  the  symplectic  structure  to  a  certain  order. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Section  2  the  symplectic  conditions  for  Runge-Kutta  methods 
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are  first  briefly  reviewed  to  fix  the  notation,  and  then  the  fixed-point  iteration-based  approximation  is  analyzed. 
Analysis  on  Newton’s  method-based  approximation  is  presented  in  Section  3.  Comparisons  among  these  approxima¬ 
tion  schemes  and  two  other  schemes  are  conducted  in  Section  4  through  various  numerical  examples  with  a  special 
focus  on  the  nonlinear  pendulum.  Finally  some  concluding  remarks  are  provided  in  Section  5. 


2  Fixed-Point  Iteration-Based  Approximation 

2.1  Symplectic  Runge-Kutta  schemes 


Consider  a  Hamiltonian  system 


Pit)  = 


(1) 


dq 

m  = 

with  the  Hamiltonian  H{p,  q),  where  {p,  q)  x  Q  for  some  integer  d  >  1,  and  Q,  the  configuration  space,  is  some 

d-dimensional  manifold.  In  this  paper  Q  =  is  assumed  for  ease  of  discussion,  but  the  extension  of  the  results  to 

/  \ 


A 


the  case  of  a  general  Q  is  straightforward.  Let  z  = 


.  Then  (1)  can  be  rewritten  as: 


\V 


A 


Kt)  =  fizit))  =  J^zH{z{t)) 


(2) 


where 


0  -Id 

Id  0 


Id  denotes  the  d-dimensional  identity  matrix,  and  Vz  stands  for  the  gradient  with  respect  to  z. 


An  s-stage  Runge-Kutta  method  to  integrate  (2)  is  as  follows  (Hairer  et  ah,  1987): 

/ 

yi  =  zo  -k  tX;Li  atjfivj),  z  =  1,  •  •  •  ,  s 

<  ,  (3) 

Zi  =  Zo  +  Tj2Ul^ifiyi) 

< 

where  t  is  the  time  step,  zq  is  the  initial  value  at  time  to,  zi  is  the  numerical  solution  at  time  to  +  t,  aij,bi  are 
appropriate  coefficients  satisfying  the  order  conditions  of  the  Runge-Kutta  method. 
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Let  'I'r  be  the  one  time-step  flow  associated  with  the  algorithm  (3),  i.e.,  zi  =  From  Sanz-Serna  (1988),  the 

transformation  'i’r  preserves  the  symplecticness  of  the  original  system  (2)  if 


^i^ij  “t”  ^j^ji  —  0,  J  —  f  j  *  *  *  5 


(4) 


Thus  if  (4)  is  satisfied,  we  have: 


OZq  OZq 


(5) 


where  stands  for  the  transpose.  The  condition  (4)  forces  the  symplectic  Runge-Kutta  method  (3)  to  be  implicit. 


To  put  (3)  in  a  more  compact  form,  denote 


f  2/1 1 


i /(2/i)l 


A 

y  = 


F(y)  = 


ysj 


yiivs)  j 


b  =  (6i,  •  •  •  ,  6s),  Aq  =  [oy  ],  and  A  =  Aq  0  hd,  where  “(8)”  denotes  the  Kronecker  (tensor)  product.  Recall  for  two 
matrices  M  =  [niij]  and  R  =  [vij],  the  Kronecker  product 


rriiiR  mi2R  ■  ■  ■ 


M®R  = 


m2iR  m22R  •  ■  ■ 


The  algorithm  (3)  can  now  be  written  as 

/ 

y  =  G(zo,  y)  =  1  O  2o  -k  rAF(y) 

(6) 

Zi  =  zo  +  Tb0l2dF(y) 
where  1  is  an  s-dimensional  column  vector  with  1  in  every  entry. 
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2.2  Approximation  based  on  fixed-point  iteration 


It  is  well-known  that  for  a  fixed  zq,  when  r  is  sufficiently  small,  there  is  a  unique  solution  y*  to  the  first  equation 
in  (6)  and  it  can  be  obtained  through  fixed-point  iteration  (Hairer  et  ah,  1987).  The  following  proposition  states  a 
similar  result;  the  key  difference  is  that  uniform  convergence  (with  respect  to  zq)  is  achieved.  As  we  shall  see,  such 
uniform  convergence  is  crucial  for  establishing  the  convergence  of  the  symplectic  structure. 

In  this  paper  ||  •  ||  will  be  used  to  denote  the  2-norm  (or  the  induced  2-norm)  of  a  vector,  a  matrix,  or  a  higher  rank 
tensor  depending  on  the  context.  For  an  open  set  fl,  its  e— neighborhood,  e),  is  defined  as 

A/’(n,  e)  =  {z  G  :  min  \\z  —  zo||  <  e}, 

Zq  G  ^ 

where  f),  denotes  the  closure  of  ft.  Denote  by  e)  the  product  of  s  copies  of  N{Q,  e), 

A7®(fl,  e)  =  Af{n,  e)  X  •  •  •  X  N{n,  e). 

Proposition  2.1  Let  Lt  C  be  a  bounded,  convex,  open  set.  Let  f  be  continuously  differentiable.  Then  for  any 
e  >  0,  there  exists  tq  >  0  dependent  on  LI  and  e  such  that,  yr  <  tq,  Vzq  G  LI, 


(1)  G(zo)’)  maps  {LI,  e)  into  itself; 

(2)  There  is  a  unique  solution  y*  to  the  first  equation  in  (6),  and  it  can  be  approximated  iteratively  via 


yN  =  G(zo,y[”-il) 
y[0l  =  1  (g) 

and 

(3)  ||yH  _y*||  <  (5"||y[0l  -y*||  with  0  <  (5  <  1,  where  5  =  tCiHAqH  and  Ci  =  max^g^(o_,)  |||^(z)||. 


(7) 


Proof.  Denote  Cq  =  maxyg_ys(f2,e)  ||F(y)||.  Let  ti  =  (note  that  ||Ao||  =  |jA||).  Then  Vr  <  n,  Vzq  G  LI, 

G(zo,-)  maps  J\f^{Ll,e)  into  itself.  Let  r2  >  0  be  such  that  r2Gi||Ao||  <  1.  Since  G(zo,-)  is  Lipschitz  continuous 
with  Lipschitz  constant  tGi||Ao||  by  the  convexity  assumption,  it  becomes  a  contraction  mapping  on  Af''^{Ll,  e)  when 
T  <  Tq  =  min{ri,r2}.  The  rest  of  the  claims  then  follows  from  the  contraction  mapping  principle  (Smart,  1974).  □ 


Remark  2.1  The  convexity  of  LI  is  assumed  only  for  using  the  mean  value  theorem  to  get  the  estimate  of  Lipschitz 
constant.  This  assumption  is  not  restrictive  since  one  can  resort  to  its  convex  hull  if  LI  is  not  convex. 
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An  explicit  but  approximate  algorithm  to  solve  (6)  is  as  follows:  for  some  >  1 


yW  =G(zo,y['=-il),  ,N 

'  y[ol=l®zo  • 


=  zo  +  T6«)/2dF(y[^l) 


From  the  implicit  function  theorem,  when  r  is  sufficiently  small,  the  solution  y*  to  the  first  equation  in  (6)  is  a 
function  of  zq,  written  as  y*(zo))  and 

^(^o)  =  [/2.d-rA^(y*(zo))]-'(l®/2d).  (9) 

ozq  ay 

Similarly  zi  in  (6),  in  (8)  (and  smooth  functions  of  them)  are  all  continuously  differentiable 

functions  of  zq.  In  the  sequel  when  we  write,  e.g.,  or  ^^F(y[^l),  we  think  of  y*  or  F(y['^l)  as  a  function  of  zq 
although  it  is  not  explicitly  written  out. 

Denote  by  the  one  time-step  flow  associated  with  the  algorithm  (8),  i.e.,  z^^^  =  The  following  lemma 

will  be  essential  for  studying  how  far  is  away  from  being  symplectic. 


Lemma  2.1  Let  Ll  C  he  bounded,  convex  and  open.  For  e  >  0,  pick  tq  as  in  the  proof  of  Proposition  2.1.  Let  f 
be  twice  continuously  differentiable  on  e).  Then  Vt  <  tq,  Vzq  G 


^  DojC!  +  CoC2N)S^+^ 
"  dzo  dzo "  -  Cf 

Ilf  (F(y™)  -  F(y))||  <  +  + 

OZq  Cl 

where  S  =  rCiHAoH, 


(10) 

(11) 


A  ^F 

Do=  max  || -  tA  — (y)]"f  1  0  72d)||  (=  max  fs|| -  tA  — (y)] 

yeA^qn,e),r<ro  Oy  yeA^qn,e),r<ro  Oy 


-1| 


Co=  max  ||F(y)||, 
yeA^®(n,e) 


^  A  ,,  9F  ,  , ,,  ,  ,,df ... 

C\=  max  h^(y)  (=  max  ), 

yeA^qn.e) "  ay  ^  zeA^(n,e) "  Sz 


(12) 


(13) 


2  =  max 

(Q,e),l<z,j'<2sd 


IQ({yij})ll,  and  Q({yi,j})  is  a  third-rank  tensor  whose  {i,j)  —  th  element  is  a  (14) 


d  9F 

vector  given  by  -^{-^)i,j{yi,j)  (here  (fp)ij  denotes  the  (i,j)  —  th  component  of  ^). 
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Proof.  See  Appendix  A.  □ 


The  main  result  of  this  section  is: 


Theorem  2.1  Let  Lt  C  be  bounded,  convex  and  open.  For  e  >  0,  pick  tq  as  in  the  proof  of  Proposition  2.1.  Let 
f  be  twice  continuously  differentiable  on  e).  Then  Vr  <  tq,  Vzq  G  LI, 

^  2\\b\\DoD,{Cf  +  CoC^il  +  N))6^+^ 

^  dzo  ’  ||Ao||C2 

,  , \\b\\Do{C!  +  CoW  +  N))6^+\^ 

^  IIAollC?  ^  ^ 


where 

A  <9F 

£>1=  max  ||/2<i  +  T&(g)/2<i^(y)[/2sd -tA  — (y)]“^(l(8)72d)||,  (16) 

yGA/ ®(n,e),r<ro  uy  uy 

and  S  and  the  other  constants  are  as  defined  in  Lemma  2.1. 


Proof.  Let  '^r  be  the  one  time-step  flow  associated  with  (6).  From  (6)  and  (8), 


A['^l(zo)  =  -  ^r(^o)  =  r&  0  J2d(F(y['^l)  -  F(y*)). 


Using  Lemma  2.1,  one  derives 


MW(^o)|I  <  r||6||ilo(Ci"  +  CoC2(l  +  A)),5^+i 
dzo  "  -  Cl 


Next  write 


i9d>i  ^(zo)  yp  II 

dzo  ’  ^  dzo  ’ 

^  AAW(^o)  ,  ^^r(^o),T  ..gAW(zo)  9vl/.(zo) 
dzo  dzo  ’  ^  dzo  dzo 

uZq  OZq  OZq 


)-J\\ 

dzo 


d^r{zo).^^T  j/d'^rjzo) 
dzo  ’  ^  dzo 


where  the  last  term  vanishes  since  T,-  is  symplectic.  The  claim  now  follows  from  (17),  ||  J||  =  1,  and 


d-^r{Zo) 

dzo 


Whd  +  rb  ®  l2d^{y*)^^\\  <  Di. 

dy  dzo 


(18) 


□ 
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Remark  2.2  Theorem  2.1  provides  a  structural  error  bound  of  in  terms  of  various  constants  specific  to  the 
problem  of  interest.  Absorbing  the  constants  and  dropping  the  second  term  in  the  right-hand  side  of  (15)(since  the 
first  term  dominates),  the  error  bound  is  simplified  to  (ci  +  for  ci,C2  >  0  and  0  <  (5  <  1.  Note  the 

connection  and  the  difference  between  this  bound  and  item  3  of  Proposition  2.1.  As  N  gets  large,  the  structural  error 
approaches  zero  and  becomes  almost  symplectic. 

3  Newton’s  Method-Based  Approximation 

Newton’s  method  is  an  alternative  to  the  fixed  point  iteration  scheme  for  solving  the  implicit  equation  in  (6).  It 
reads 

yW  =  G(zo,y'"-'')  =  y'"-"'  -  [hsd  -  TA^(y["-il)]-i(y["-il  -  1  0  zq  -  rAF(y["-il)).  (19) 

Typically  convergence  conditions  for  Newton’s  method  include  that  the  Jacobian  is  invertible  at  the  solution  point 
and  that  the  initial  condition  is  close  enough  to  the  solution  (Schwarz,  1989).  Such  conditions  often  cannot  be  verified 
directly.  For  the  special  case  (6),  however,  Proposition  3.1  shows  that  when  taking  the  natural  candidate  for  y[°l, 
the  convergence  is  guaranteed  if  r  <  tq,  where  tq  can  be  determined  explicitly. 

Proposition  3.1  Let  Lt  C  be  a  bounded,  convex,  open  set.  Let  f  be  three  times  continuously  differentiable.  Then 
for  any  e  >  0,  there  exists  tq  >  0  dependent  on  LI  and  e  such  that,  yr  <  tq,  Vzq  G  LI, 

(1)  G(zo)’)  maps  {LI,  e)  into  itself; 

(2)  There  is  a  unique  solution  y*  to  the  first  equation  in  (6),  and  it  can  be  approximated  iteratively  via 

/ 

yN  =  G(zo,y["-^l) 

<  ;  (20) 

y[0l  =  1  (g)  zo 

and 

(3)  ||y["'  -yll  <  -yiP",  where  K>0  and  A||y*  -y[ol||  <  1. 


Proof.  Through  algebraic  manipulations,  G(2:o)y)  can  be  rewritten  as 


Pick  Ti  >  0  such  that  l2sd  —  is  invertible  Vr  <  ti,  Vy  G  A/’^(n,  e).  Let 

Eq=  max  ||[/2sd -TA^(y)]"i||, 

(0,e),T<Ti  dy 

„ll|^(y)(i«>^o-y)  +  F(y)||, 
y&M‘(Q.,i),zo&n  ay 


(22) 

(23) 


and  let  T2  >  0  be  such  that  T2E0E1 1|  Ao||  <  1.  Then  it  can  be  verified  that  if  t  <  min{Ti,  r2},  6(2:0,  •)  maps  e) 

into  itself. 


The  next  goal  is  to  establish  that  G(zo,  •)  is  a  contraction  mapping.  This  can  be  done  by  evaluating  To  properly 
handle  the  third-rank  tensor  involved,  for  rj  G  one  calculates 

^2-p 

-^(^o.y)  •  V  =  -'rH(y)A(^(y)  •  77)H(y)[y  -  1  0  zq  -  'rAF(y)],  (24) 

where  denotes  the  action  of  a  second-rank  or  third-rank  tensor  on  a  vector,  and 

H(y)  =  [/2.d-TA^(y)]-i.  (25) 

Denote 


E2=  niax  ||^(y)||, 
yeM‘iQ,e)  dy^ 

E3  =  max  ||y 

(n,e),2oG^,r<Ti 


l(8)2;o-rAF(y)||, 


(26) 

(27) 


and  pick  T3  >  0  such  that  T3EqE2E3\\Ao\\  <  1.  Then  when  r  <  min{ri,  r2,  T3},  6(20,  •)  is  a  contraction  mapping 
and  hence  (20)  converges  to  a  (unique)  fixed  point,  which  is  the  solution  to  the  first  equation  in  (6). 


Since  ^{zq,  y*)  =  0,  the  convergence  rate  of  (20)  is  quadratic,  as  is  standard  for  Newton’s  method  (Schwarz,  1989): 


IlyN  -y*||  <  if|jy["-l]  _y*||2  <  "  1  1 1  y  M  -  y  *  f "  , 


(28) 


where 


.d^G 


max 


(^o,y)||- 


(29) 


yGA/^^(n,e),2o£f2,'F<Ti  dy^ 

It’s  easy  to  see  that  ^^(20, y)  contains  a  factor  of  t.  On  the  other  hand,  ||yt°l  —  y*||  <  tCoHAoH,  where  Cq  is  as 
defined  in  Lemma  2.1.  Therefore  there  exists  r4  >  0  such  that  when  t  <  A|jy*  —  y[°]||  <  1.  Finally  tq  in  the 
statement  of  the  proposition  can  be  chosen  to  be  tq  =  minjri,  r2,  T3,  T4}.  □ 
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Analogous  to  (8),  an  approximation  scheme  for  solving  (6)  can  be  constructed  based  on  Newton’s  method:  for  some 

N>1, 


.W  =G(0o,y''=-''),  ,N 


.[0]  =  1 . 


=  zo  +  T6«)/2dF(y[^l) 


i  Zo 


(30) 


Denote  by  the  one  time-step  flow  associated  with  the  algorithm  (30).  The  following  two  lemmas  will  be  used 
in  the  proof  of  Theorem  3.1. 


Lemma  3.1  Let  LI  C  be  bounded,  convex  and  open.  For  e  >  0,  pick  tq  as  in  the  proof  of  Proposition  3.1.  Let 
f  be  three  times  continuously  differentiable  on  J\f{fl,e).  Define  H(-)  as  in  (25),  and  J(y)  =  H(y)A^(y).  Then 
Vt  <  To,  Vzo  G  LI, 


II  „  II  <  Gy  —  v^(l -I- - - ), 

ozo  1  -  7o 

|AH(y[^l)||<Cff 

OZo  E3 

||_^j(y[Af])||  <  =  II^o||(Gi7o  -I-  EoE2E3)Cy 

dzo  ~  E3 


A 


(31) 

(32) 

(33) 


where  70  =  T3EqE2E3\\A3\\;  Cf  is  as  defined  in  (13);  Ei,  E2  are  as  defined  in  (23),  (26);  and  Eg  and  E3  are  as 
defined  in  (22),  (21)  with  ti  replaced  by  tq. 


Proof.  See  Appendix  B.  □ 


Lemma  3.2  Let  LI  C  be  bounded,  convex  and  open.  For  e  >  0,  pick  tq  as  in  the  proof  of  Proposition  3.1.  Let  f 
be  three  times  continuously  differentiable  on  Af(Ll,  e).  Then  Vt  <  tq,  Vzq  G  LI, 


dzn  dzf) 


(34) 

(35) 


where  6  =  tGo||Ao||A'  <  1,  Dy  =  ^{Cj  +  GiGTr||Ao||  -I-  ^G2Dg||Ao||),  Cj  and  Ch  are  as  defined  in  Lemma  3.1, 
and  Cl,  C2,  Dq  and  K  are  as  defined  in  (13),  (1)),  (12)  and  (29),  respectively. 


A 


Proof.  See  Appendix  C.  □ 
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Following  the  arguments  as  in  the  proof  of  Theorem  2.1  and  using  Lemma  3.2,  we  can  show: 


Theorem  3.1  Let  Lt  C  be  bounded,  convex  and  open.  For  e  >  0,  pick  tq  as  in  the  proof  of  Proposition  3.1. 
Let  f  be  three  times  continuously  differentiable  on  Af{Ll,e).  Let  be  the  one  time-step  flow  associated  with  (30). 


Then  Vr  <  tq,  V^q  €  LI, 


)  -  J\\  <2TDi\\b\\{CiDy6'^''~" 


C2D0 

K 


5^  )  +  {T\\b\\{C^Dy5^ 


C2D0 

K 


<5^"))"(36) 


where  Di  is  as  defined  in  (16),  and  6  and  the  other  constants  are  as  defined  in  Lemma  3.2. 


4  Numerical  Examples  and  Discussion 

The  performances  of  approximation  schemes  (8)  and  (30)  on  symplectic  structure  conservation  have  been  character¬ 
ized  in  Theorem  2.1  and  Theorem  3.1,  respectively.  Under  suitable  conditions  and  with  proper  choices  for  the  step 
size  and  the  initial  iteration  value  both  schemes  uniformly  (with  respect  to  zq)  converge,  and  the  convergence 
rate  of  symplectic  structure  for  either  scheme  is  closely  connected  to  the  corresponding  rate  for  the  solution  conver¬ 
gence  (i.e.,  —  y*||)-  In  this  section  the  implications  of  these  results  for  implementing  symplectic  Runge-Kutta 

schemes  are  explored  through  a  variety  of  numerical  examples. 

Important  factors  in  choosing  a  Runge-Kutta  scheme  for  Hamiltonian  systems  include  the  numerical  accuracy, 
the  structural  preservation  performance  (symplecticness)  and  the  computational  cost.  Since  the  issue  of  numerical 
accuracy  is  not  the  focus  of  this  paper,  the  discussion  will  be  centered  around  the  interplay  between  the  symplecticness 
and  the  computational  complexity.  For  illustrational  purposes,  the  methods  listed  in  Table  1  will  be  compared  in  the 
numerical  problems.  For  a  definition  of  pseudo-symplecticness  order,  we  refer  to  Aubry  and  Chartier  (1998).  The 
mid-point  rule  and  the  Gauss  method  are  implicit,  and  both  fixed-point  iteration  and  Newton’s  method  will  be  used 
to  solve  the  implicit  equations.  Table  2  lists  the  test  problems.  Some  of  these  problems  were  also  used  in  Aubry  and 
Chartier  (1998).  The  computation  was  done  in  Matlab  on  a  Dell  laptop  Inspiron  4150. 

4.1  The  nonlinear  pendulum  problem 

An  essential  property  of  a  symplectic  map  is  area-preservation.  The  ellipse  shown  in  Fig.  1,  with  semi-major  axis  =1.8 
and  semi-minor  axis  =  1.2,  represents  the  set  of  initial  conditions  for  the  nonlinear  pendulum  problem.  Numerical 
solutions  after  one  time-step  under  different  methods  are  compared  with  the  exact  solution  in  Fig.  2,  where  the  time 
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Table  1 


Runge-Kutta  methods  used  in  numerical  examples. 


Notation 

Method 

Order 

Pseudo-symp.  order 

s 

MidPoint 

Mid-point  rule 

2 

symplectic 

1 

Gauss4 

Gauss  method  (Hairer  et  al.,  2002) 

4 

symplectic 

2 

PS63 

Pseudo-symp.  method  (Aubry  and  Chartier,  1998) 

3 

6 

5 

RK4 

Classical  Runge-Kutta 

4 

4 

4 

Table  2 

Test  problems  used  in  the  numerical  study. 


Problem 

Hamiltonian  H {p,  q) 

Step  size  r 

Initial  Conditions 

Nonlinear  pendulum 

4  -  cos(<j) 

1.6,  0.8,  0.2 

See  the  text 

Linear  pendulum 

0.5 

(2,2)^ 

Kepler  problem 

0.1 

(0,2,  0.4,0)"^ 

Bead  on  a  wire 

1  TTf^\  ..rUU  TTf^\ 

=  0.1(g((7  - 

1 

(0.49,  0)"^ 

2(l  +  U'(qy-‘)  ^  ^  1*'-'  V  (q) 

2)f  +  0.008(7® 

6 

Galactic  dynamics 

5  (Pl  +P2  +P3)  +  jiPl(l2  — P2(7l)-|-ln(l  -1  ^  -f 

fi-  +  ^),  with  a=|,6=l,c=| 

0.2 

(0,1.689,0.2,2.5,0,0)'^ 

step  r  =  1.6,  and  the  implicit  equations  in  MidPoint  and  Gauss4  were  solved  using  Newton’s  method  up  to  machine 
accuracy.  As  one  can  see,  the  (exact)  final  configuration  is  distorted  from  the  initial  elliptical  curve.  By  symplecticity 
of  the  exact  flow,  the  area  enclosed  by  the  exact  solutions  at  t  =  1.6  is  equal  to  that  enclosed  by  the  initial  curve. 
Among  the  numerical  solutions,  Gauss4  has  the  best  performance  in  terms  of  accuracy  and  area-preservation  since 
it  completely  overlaps  the  exact  solution.  The  solution  of  MidPoint  is  noticeably  different  from  that  of  the  exact  one 
because  it  is  of  the  second  order.  The  area-preserving  performance  of  MidPoint  cannot  be  easily  told  from  the  figure 
(theoretically  it  should  be  as  good  as  that  of  Gauss4).  Under  PS63  it  can  be  seen  that  the  area  has  shrunk  a  little 
bit,  while  RK4  delivers  the  worst  performance  in  area  preservation. 

To  provide  a  quantitative  measure  of  area  preservation,  we  have  picked  10^  points  on  the  ellipse  (as  initial  conditions). 
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The  initial  area  Ao  at  time  t  =  0  is  approximated  by  the  sum  of  areas  of  10^  triangles  formed  by  the  picked  points 
and  the  origin  (see  Fig.  1  for  illustration,  where  8  points  are  used).  The  final  area  at  <  =  1.6  is  calculated  similarly, 
using  the  current  10^  solution  points.  The  (normalized)  area  error  is  then  defined  as  . 


One  goal  of  this  paper  is  to  provide  insight  into  the  choice  of  fixed-point  iteration  versus  Newton’s  method.  From 
Theorem  2.1  and  Theorem  3.1,  Newton’s  method  enjoys  much  faster  structural  convergence  than  the  fixed-point 
iteration  in  terms  of  the  number  of  iterations.  This  is  verified  in  Fig.  3  and  Fig.  4.  Fig.  3  shows  the  decrease  of  area 
error  with  the  number  of  fixed-point  iterations,  where  the  underlying  algorithm  used  was  MidPoint.  In  the  figure, 
the  bound  from  Theorem  2.1  is  also  plotted.  Note  the  similar  trend  in  both  curves,  in  particular,  their  consistent 
convergence  rates.  In  Fig.  4,  the  area  error  stops  decreasing  after  4  iterations,  and  it  stays  around  0.5  x  10“®.  This 
is  not  due  to  the  limitation  of  machine  accuracy.  It  is  considered  to  arise  from  approximating  the  area  based  on  10^ 
points.  Indeed,  we  have  observed  that  the  error  stops  decreasing  around  10“®  if  10^  points  are  used  in  computing 
the  area. 


Despite  the  faster  convergence,  Newton’s  method  takes  longer  time  in  each  iteration  than  the  fixed-point  iteration. 
This  brings  up  the  issue  whether  the  aforementioned  advantage  is  still  an  advantage  when  actual  computation  time  is 
considered.  In  terms  of  N,  the  computation  times  of  the  two  methods  can  be  approximately  expressed  as  Tq  +  NTf, 
Tg  -I-  NTi,  respectively.  Here  Tg  and  Tf  represent  the  computation  overhead  and  the  computation  cost  per  iteration 
for  the  fixed-point  scheme,  respectively,  and  Tg  and  Tf  represent  the  counterparts  for  Newton’s  method.  The  actual 
computation  times  taken  by  the  two  methods  are  plotted  in  Fig.  5,  both  displaying  a  linearly  increasing  trend.  As 
N  gets  large,  the  ratio  of  their  computation  costs  approaches  a  constant  Considering  their  convergence  rates, 
one  can  conclude  that  Newton’s  method  is  more  time-efficient  when  very  low  structural  error  is  needed. 

Two  other  step  sizes  r  =  0.8,  r  =  0.2  are  used  to  integrate  the  nonlinear  pendulum  equation  while  the  final  time  is 
kept  the  same,  i.e.,  t  =  1.6.  Therefore  the  time  steps  for  these  step  sizes  are  2  and  8,  respectively.  Fig.  6  shows  the 
work-precision  diagrams  of  the  fixed-point  iteration  scheme  for  the  three  different  step  sizes.  It  can  be  seen  that  for 
the  same  amount  of  CPU  time,  with  t  =  0.8,  the  area  error  is  smaller  than  that  with  r  =  1.6  or  with  r  =  0.2.  It 
can  be  explained  as  follows:  when  r  is  relatively  big,  the  convergence  rate  is  slow;  while  when  r  is  relatively  small, 
it  requires  many  time  steps  which,  to  keep  the  total  CPU  time  the  same,  leads  to  a  small  number  N  of  iterations  at 
each  time  step.  Therefore  to  maximize  the  computational  efficiency  (defined  as  the  level  of  structural  preservation 
per  CPU  time  unit),  one  needs  to  seek  a  moderate  step  size.  Fig.  7  shows  the  work-precision  diagrams  of  Newton’s 
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method-based  approximation  under  different  step  sizes.  For  this  particular  problem,  even  with  r  =  1.6,  at  most  4 
iterations  would  bring  the  area  error  down  to  the  order  of  10“®  (the  achievable  limit  as  explained  earlier),  and  there 
is  not  much  to  gain  by  using  smaller  r  in  the  sense  of  computational  efficiency  defined  above. 

Fig.  8  through  Fig.  10  compare  the  work-precision  diagrams  of  Gauss4/FixedPt  (solving  Gauss4  with  fixed-point 
iteration),  Gauss4/Newton  (solving  Gauss4  with  Newton’s  method),  PS63  and  RK4  for  different  step  sizes.  PS63 
always  beats  RK4  at  a  slight  cost  of  computational  time.  For  same  amount  of  GPU  time,  PS63  also  leads  to  smaller 
area  error  than  Gauss4/FixedPt  and  Gauss4/Newton.  However,  while  the  structural  error  under  Gauss4/FixedPt  or 
Gauss4/Newton  approaches  zero  with  increasing  GPU  time,  the  error  of  PS63  can  be  large  when  r  is  relatively  big 
(Fig.  8  and  Fig.  9).  Finally,  it  can  be  seen  that  corresponding  to  relatively  large  error,  Gauss4/FixedPt  needs  less  GPU 
time  than  Gauss4/Newton;  but  for  very  small  error,  Gauss4/Newton  requires  less  GPU  time  than  Gauss4/FixedPt. 

From  (28)  and  the  proof  of  Lemma  3.2,  a  better  choice  of  (i.e.,  smaller  —  y*||  with  y[°l  smoothly  dependent 
on  zq)  leads  to  faster  convergence  of  the  symplectic  structure.  A  hybrid  approximation  scheme  is  motivated  by  this 
observation:  first  use  iGzg  as  the  initial  guess  and  run  the  fixed-point  iteration  Ni  times,  then  use  as  the  initial 
value  and  run  Newton’s  method  for  N2  iterations.  The  idea  is  to  use  relatively  cheap  computation  of  the  fixed-point 
algorithm  to  get  a  better  initial  estimate  for  Newton’s  method.  Fig.  11  shows  the  work-precision  comparison  of 
this  hybrid  scheme  (with  Ni  =  1)  with  the  plain  Newton’s  method,  where  both  cases  of  r  =  1.6  and  t  =  0.8  are 
displayed.  From  the  figure  it  can  be  seen  that  the  hybrid  scheme  offers  faster  convergence  rate  with  a  slight  increase 
of  computational  cost.  Again  when  looking  at  the  figure,  one  should  keep  in  mind  that  0.5  x  10“®  is  the  achievable 
area-error  limit  as  a  result  of  our  area-approximation  method,  and  hence  he  or  she  should  not  be  confused  by  the 
somewhat  misleading  slopes  of  the  last  segments  of  the  curves. 

4-2  Other  problems 

Fig.  12  shows  the  trajectories  of  the  linear  pendulum  in  the  phase  space  under  Gauss4/FixedPt  and  Gauss4/Newton 
for  5  X  10^  time  steps  (the  data  were  down-sampled  by  20  to  reduce  the  file  size).  For  the  fixed-point  iteration 
method,  the  energy  decays  to  zero  if  the  iteration  number  N  =  3.  The  energy  decay  rate  is  significantly  reduced 
when  N  is  increased  to  5,  and  with  N  =  8,  the  trajectory  almost  stays  on  the  circular  orbit.  Newton’s  method,  on 
the  other  hand,  gives  rise  to  the  exact  solution  (up  to  the  machine  precision)  in  one  iteration  since  the  system  is 
linear. 
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The  numerical  solutions  of  the  Kepler  problem  (qi  and  q2  components)  are  plotted  in  Fig.  13  (after  down-sampling  by 
20).  For  comparison,  the  exact  orbit  is  also  shown.  The  total  number  of  time  steps  is  2  x  10^.  It  can  be  observed  that 
when  N  =  6,  the  solution  with  Gauss4/FixedPt  follows  a  precession  motion  of  the  elliptical  orbit.  This  is  also  true  for 
Gauss4/Newton  {N  =  2).  Such  “precession”  effect  is  typical  when  integrating  the  Kepler  problem  with  a  symplectic 
scheme  (Hairer  et  ah,  2002).  Note  that  PS63  also  demonstrates  the  similar  behavior  with  a  slower  precession  rate. 
For  Gauss4/FixedPt  with  fV  =  2,4  and  RK4,  the  solutions  distort  the  ellipse.  The  angular  momentum  is  also  a 
conserved  quantity  for  the  Kepler  problem.  Fig.  14  shows  the  angular  momentum  error  under  different  schemes. 
Listed  in  Table  3  is  the  GPU  time  used  in  the  computation. 

Fig.  15  and  Fig.  16  show  the  evolution  of  error  in  the  Hamiltonian  for  the  bead-on-a-wire  problem  and  the  galactic 
dynamics  problem.  Table  4  and  Table  5  list  the  GPU  time  used  by  different  algorithms. 

Table  3 

CPU  time  used  in  solving  the  Kepler  problem. 


Method 

Gauss4/FixedPt 

Gauss4 /N  ewton 

PS63 

RK4 

N 

2 

4 

6 

8 

N 

2 

3 

Time  (sec.) 

49.9 

70.6 

89.3 

108.2 

73.0 

91.0 

55.0 

44.7 

Table  4 

CPU  time  used  in  solving  the  bead-on-a-wire  problem. 


Method 

Gauss4/FixedPt 

Gauss4 /N  ewton 

PS63 

RK4 

N 

3 

5 

7 

N 

1 

2 

3 

Time  (sec.) 

160.4 

197.0 

240.0 

158.3 

191.0 

220.4 

146.4 

129.2 

Table  5 

CPU  time  used  in  solving  the  galactic  dynamics  problem. 


Method 

Gauss4 /FixedPt 

Gauss4/Newton 

PS63 

RK4 

N 

2 

4 

N 

1 

2 

Time  (sec.) 

283.9 

360.0 

311.9 

376.4 

318.2 

264.2 

15 


5  Conclusions 


Symplectic  Runge-Kutta  schemes  for  the  integration  of  general  Hamiltonian  systems  are  implicit.  When  approxi¬ 
mation  methods  are  used  to  solve  the  implicit  equations,  the  resulting  integration  schemes  do  not  fully  preserve 
the  symplectic  structure  of  the  original  systems.  It  is  thus  of  interest  to  understand  the  structural  error  incurred 
by  the  approximation  schemes.  In  this  paper  approximations  based  on  two  common  iterative  methods  for  solving 
implicit  equations,  fixed-point  iteration  and  Newton’s  method,  were  analyzed  and  the  corresponding  error  bounds 
established.  Under  proper  conditions,  these  schemes  become  almost  symplectic  as  the  iteration  number  N  gets  large. 
Although  the  results  show  that  the  structural  convergence  of  either  scheme  is  closely  related  to  its  numerical  con¬ 
vergence,  the  former  (essentially  —  ^;7ll)  does  not  follow  merely  from  the  latter  (||y[^^  —  y*||);  instead  it 

is  a  consequence  of  the  uniform  convergence  of  the  iterative  schemes  with  respect  to  the  initial  condition  zq,  the 
particular  choices  of  the  initial  iteration  values,  and  the  smoothness  of  the  mappings  G(-,  •)  and  G(-,  •). 

The  theoretical  results  can  be  used  in  selecting  an  appropriate  approximation  scheme  when  integrating  a  specific 
problem.  The  emphasis  here  is  the  trade-off  between  the  computational  cost  and  the  structural  preservation  perfor¬ 
mance  although  the  numerical  accuracy  (the  order  of  a  method)  also  plays  an  important  role  in  implementation.  The 
faster  convergence  rate  of  Newton’s  method-based  scheme  makes  it  more  favorable  than  the  fixed-point  iteration- 
based  scheme,  especially  when  very  small  structural  error  is  required.  This  was  verified  in  the  numerical  tests. 

The  effect  of  the  step  size  on  the  computational  efficiency  was  studied  in  the  numerical  experiments.  We  also  note 
that  the  arguments  in  the  proofs  of  Proposition  2.1  and  Proposition  3.1  may  be  used  to  find  the  step  size  tq  (below 
which  the  scheme  is  convergent)  for  the  specific  problem  of  interest.  For  stiff  problems,  tq  will  be  very  small  for 
the  fixed-point  algorithm  and  Newton’s  method  is  generally  more  efficient.  After  observing  that  a  better  initial 
guess  would  speed  up  the  convergence  rate  of  Newton’s  method,  a  hybrid  scheme  (running  one  or  several  fixed-point 
iterations  to  obtain  initial  values  for  Newton’s  method)  was  proposed  and  explored.  Simulation  suggested  that  the 
hybrid  scheme  has  a  potential  to  out-perform  the  plain  Newton’s  method. 

The  almost  symplectic  schemes  were  also  compared  against  a  pseudo-symplectic  method  and  a  non-symplectic 
method.  It  is  of  no  surprise  that  the  non-symplectic  method  delivers  the  poorest  performance  in  area-conservation 
and  energy-conservation.  For  methods  of  comparable  orders  of  accuracy,  the  pseudo-symplectic  one  delivers  slightly 
better  structural  preserving  performance  than  an  approximation-based  symplectic  scheme  if  the  latter  spends  the 
same  amount  of  CPU  time.  However,  with  increased  CPU  time  (which  is  still  comparable  to  the  CPU  time  used 
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by  the  pseudo-symplectic  one)  the  approximation  scheme  has  the  potential  to  reach  very  low  structural  error  and 
becomes  almost  symplectic.  On  the  other  hand,  as  admitted  in  Aubry  and  Chartier  (1998),  the  design  of  a  pseudo- 
symplectic  method  (in  particular,  of  order  p  and  of  pseudo-symplecticness  order  2p  (Aubry  and  Chartier,  1998)) 
beyond  order  (3,6)  is  very  complicated.  This  will  hinder  the  use  of  pseudo-symplectic  methods  in  very  long  time 
simulation  of  Hamiltonian  systems. 
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A  Proof  of  Lemma  2.1 


Proof.  From  (6)  and  (8), 


yW  _  y*  =  rA(F(y['^-il)  -  F(y*)). 


Taking  derivative  on  both  sides  of  (A.l)  with  respect  to  zq  and  re-arranging  terms,  one  gets 


ozo  ay  dy  dzo 


(A.l) 


(A.2) 


Eq.  (A.2)  implies 


I - 1^11  <  .Moll  11^  -  II  +  y||.4„||  II -  A(y-)||  llff I 

ozq  azo  '^y  azo  '^y  ^y  <9^0 

<rCi||Ao||  11^^^  -  1^11  +rilo||7lo||  |||^(y'^-'')  -  ^(y*)l|. 
ozq  ozq  oy  ay 


(A.3) 


By  the  mean  value  theorem,  the  (t,  j)— th  component  of  fp(y^'^  ^^)  —  ^^(y*)  can  be  expressed  as 


dy' 


f) 

•  (y'^”^'  -y*)  for  some  yi^  G  A7^(U,e), 


which  leads  to 

||^(y'^-^')-|^(y*)||<C2||yfo-il-y*|| 

<  C2(5^"ly[°'  -  y*||  (from  Proposition  2.1) 

<  'rC'oC'2||Ao||(5^"^||  (since  y*  -  yf°l  =  rAF(y*)).  (A.4) 


Plugging  (A.4)  into  (A.3)  and  performing  recursions,  one  has 

II  5y*ii  .  dy*  ,  ,  CoC2DoN5^+^ 

"  dzo  dzQ "  -  "  dzo  dzQ "  ^  Cl 
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Eq.  (10)  is  then  proved  by  noting 


II 


dzQ 


(9F  9v* 

\\rA—{y*)^J<TC^Do\\A4. 


To  show  (11),  write 


A(F(y[^l)-F(y*)) 

OZq 


OF 

dy 


ay[W] 

dzo 


ozo  dy  dy  dzo 


and  then  use  (10)  and  (A. 4).  □ 


(A.6) 


B  Proof  of  Lemma  3.1 


Proof.  Differentiating  both  sides  of  (19)  with  respect  to  zq  leads  to 


5y[^i  ac,  ,  ac, 

—  =  i^(^o,y‘  ‘)^^ - +  7^(^o,y‘  ')• 

9zo  dy  OZq  OZq 


(B.l) 


From  |||p(2o,y^^  ^')||  <  Ao||  (recall  (24))  and  |||£(^:o,y^^  ^')||  =  ||H(y[^  il)l(8)72<i||  <  ^/sAo,  one  gets 


dzQ 


1^11  <rA2E2  As  Poll  11^^^ - ll  +  yiAo. 

dzo  dzo 


Since 


dzo 


A 


—  y  s, 


1%'^'  II  ^  /Ty.  JV  ,  Ao(l  -7^) 


dzn 


<  + 


1-7 


), 


where  7  =  tEqE2E3\\Ao\\.  Eq.  (31)  then  follows  from  0  <  7  <  70  <  1. 


Eq.  (32)  can  be  shown  by  writing 

^H(yl»l)  =  Aff  =  rH(yl»l)A^(yl»l)H(yl»l)  pi 

azo  dy  dzo  dy^  dzo 


and  then  using  (31). 


Finally  to  show  (33),  note  that 


#-J(yl»l)  =  AH(yM)A  fl(yl»l) 

dzQ  dzo  dy 


H(yI~l)A|lI(yl»l)-'^'"' 


■ay2 


azo 


and  then  use  (31)  and  (32).  □ 
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C  Proof  of  Lemma  3.2 


Proof.  From  (19)  and  y*  =  G(zo,y*)  ,  one  can  derive 


yim  _  y*  =  _^j(y[^-i])(y[w-i]  _  y*)  +  T H ( y ^ ) A ( F (y ^ )  -  F(y*)).  (C.l) 


Taking  derivative  on  both  sides  of  (C.l)  with  respect  to  zq  and  using  (A. 6),  it  can  be  shown  that 


dyW 

dzo 


dy* 

dzo 


-T^J(y'"-‘l)(yl"-'l  -  y*)  +  t  AH(yl»->l)A(F(yt'*-'l)  -  F(y-)) 


By  Lemma  3.1  and  the  mean  value  theorem,  Eq.  (C.2)  implies  that 

11%^^  -  <  r{Cj  +  C,Ch\\Ao\\  +  J-C2Dl\\Ao\\)\\y^^-^^  -  y*||. 

OZo  dZo  y'S 

Eq.  (34)  then  follows  from  Proposition  3.1.  Eq.  (35)  is  obtained  by  making  use  of  (A. 6)  and  (34).  □ 


q 


(C.2) 


Fig.  1.  Initial  conditions  for  the  nonlinear  pendulum  problem  and  the  schematic  of  approximating  the  enclosed  area  with  a 
finite  number  of  triangles. 
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Fig.  2.  Comparison  of  numerical  solutions  with  the  exact  one  at  t  =  1.6  (r  =  1.6)  for  the  nonlinear  pendulum  problem. 


Fig.  3.  Decrease  of  the  area  error  vs  the  number  N  of  iterations  with  fixed-point  iteration  computed  for  the  nonlinear  pendulum 
problem.  MidPoint  is  used  with  r  =  1.6  and  the  number  of  time  steps  is  one. 
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Fig.  4.  Decrease  of  the  area  error  vs  the  number  N  of  iterations  with  Newton’s  method  computed  for  the  nonlinear  pendulum 
problem..  MidPoint  is  used  with  r  =  1.6  and  the  number  of  time  steps  is  one. 


Fig.  5.  Comparison  of  the  computation  time  (for  one  time-step)  vs  the  number  N  of  iterations  for  fixed-point  iteration  and 
Newton’s  method.  The  nonlinear  pendulum  problem  is  computed  and  MidPoint  used  with  r  =  1.6. 
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Fig.  6.  Work-precision  diagrams  for  the  nonlinear  pendulum  problem  under  the  fixed-point  iteration  scheme  with  different 
step  sizes.  Final  time  t  =  1.6  fixed.  Underlying  algorithm:  MidPoint. 


Fig.  7.  Work-precision  diagrams  for  the  nonlinear  pendulum  problem  under  Newton’s  method-based  scheme  with  different 
step  sizes.  Final  time  t  =  1.6  fixed.  Underlying  algorithm:  MidPoint. 
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Fig.  8.  Comparison  of  work-precision  diagrams  for  the  nonlinear  pendulum  problem  under  different  schemes  (t  =  1.6).  Final 
time  t  =  1.6. 
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Fig.  9.  Comparison  of  work-precision  diagrams  for  the  nonlinear  pendulum  problem  under  different  schemes  (t  =  0.8).  Final 


time  t  =  1.6. 
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Fig.  10.  Comparison  of  work-precision  diagrams  for  the  nonlinear  pendulum  problems  under  different  schemes  (r  =  0.2).  Final 
time  t  =  1.6. 
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Fig.  11.  Comparison  of  work-precision  diagrams  for  the  nonlinear  pendulum  problem  under  the  hybrid  scheme  and  Newton’s 
method,  where  the  underlying  algorithm  is  Gauss4.  Gauss4/Hybrid:  run  fixed  point  iteration  once  and  then  run  Newton’s 
method,  (a)  r  =  1.6;  (b)  t  =  0.8.  Final  time  t  —  1.6  for  both  (a)  and  (b). 
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Fig.  12.  Trajectories  of  the  linear  pendulum  in  the  phase  space  under  Gauss4/FixedPt  and  Gauss4/Newton. 
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Fig.  13.  Exact  and  numerical  solutions  of  the  Kepler  problem.  The  underlying  algorithm  for  FixedPt  and  Newton  was  Gauss4. 
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Hamiltonian  error 


Fig.  14.  Comparison  of  the  angular  momentum  error  for  the  Kepler  problem. 
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Fig.  15.  Comparison  of  the  Hamiltonian  error  for  the  bead-on-a-wire  problem. 
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Fig.  16.  Comparison  of  the  Hamiltonian  error  for  the  galactic  dynamics  problem. 
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